package fem2.material;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import fem2.Element;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Model;
import fem2.Node;
import fem2.Patch;
import fem2.PatchManager;
import fem2.element.StructuralElement;
import fem2.enu.StressType;

/**
 * A patch stress recovery formulated as a host material model using Zhu &
 * Zienkiewicz estimator
 * 
 * @author hbui
 * 
 */
public class SPRrecoveryMM extends StressRecoveryMM {

	private PatchManager pm;

	public SPRrecoveryMM(Model m, MaterialModel mm) {
		this(new PatchManager(m.getMesh()), mm);
	}

	public SPRrecoveryMM(PatchManager pm, MaterialModel mm) {
		super(mm);
		this.pm = pm;
	}

	public PatchManager getPatchManager() {
		return this.pm;
	}

	@Override
	public IVector getStress(Element e, int k, StressType stressType) {
		switch (stressType) {
		case RECOVER_STRESS_ZZ:
			return computeRecoverStress(e, e.getIntegrationPoint(k));
		default:
			return getMaterialModel().getStress(e, k, stressType);
		}
	}

	/**
	 * compute recover stress at node
	 * 
	 * @param n
	 * @return
	 */
	public IVector computeRecoverStress(Node n) {
		ArrayList<Patch> patches = pm.getPatch(n);
		IVector rStress = null;

		for (int i = 0; i < patches.size(); i++) {
			IVector tmp = computeRecoverStress(n, patches.get(i));
			if (rStress == null) {
				rStress = tmp;
			} else {
				BLAM.add(1.0, tmp, rStress);
			}
		}

		BLAM.scale(1.0 / patches.size(), rStress);

		return rStress;
	}

	/**
	 * compute contribution of a patch to recover stress at a node
	 * 
	 * @param n
	 * @param ntrs
	 * @param p
	 */
	private IVector computeRecoverStress(Node n, Patch p) {
		int dim = n.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		ArrayList<IVector> stressList = new ArrayList<IVector>();
		IVector rStress = new ArrayVector(ntrs);
		/*
		 * gather all stress at Gauss point around the node
		 */
		for (int i = 0; i < p.countMeshPart(); i++) {
			MeshPart mp = p.getMeshPart(i);
			StructuralElement se = (StructuralElement) mp.getElement();
			for (int j = 0; j < se.getNumIntegrationPoints(); j++) {
				stressList.add(se.getStress(j));
			}
		}
		/*
		 * collect each stress component to form coefficients and solve for
		 * respective recovered stress
		 */
		double[] values = new double[stressList.size()];
		for (int i = 0; i < ntrs; i++) {
			for (int j = 0; j < stressList.size(); j++) {
				values[j] = stressList.get(j).get(i);
			}
			IVector coeff = p.solveForCoefficients(values);
			IVector P = p.computeInterpolationVector(n.getPosition());

			// Debugger.watch("coeff = ", coeff);
			// Debugger.watch("P = ", P);

			rStress.set(i, BLAM.dot(P, coeff));
		}
		return rStress;
	}
}
